
#define SOL_ALL_SAFETIES_ON 1
#include <sol/sol.hpp>
#include "matrix.hpp"
#include "eigen_sparse.h"
#include "eigen_sparse_spectra.h"
#include <Eigen/Sparse>
//#include <Eigen/SparseCore>
#include <iostream>
#include <list>
#include <cmath>
#include <complex>

typedef Eigen::SparseMatrix<double> SpMat;
typedef Eigen::SparseMatrix<std::complex<double>> SpCMat;
typedef Eigen::VectorXcd EigenVectorXcd;
namespace sol {
        //template <>
        //struct is_automagical<EigenMatrix> : std::false_type {};
        //template <>
        //struct is_automagical<EigenVectorXcd> : std::false_type {};
		template <>
		struct is_to_stringable<SpMat> : std::true_type {};
}

template< class T >
std::ostream& to_string1(std::ostream& os, const Eigen::SparseMatrix<T>& right) {

	os << "sp_matrix" <<std::endl;
	if (right.size() == 0) {
		os << "{ sp_matrix! empty }";
		return os;
	}

	os << "{ sp_matrix " << right.size() << " | ";

	for (int k=0; k<right.outerSize(); ++k)
	  for (SpMat::InnerIterator it(right,k); it; ++it)
	  {

	    os << "[" << it.row();   // row index
	    os << "," << it.col();   // col index (here it is equal to k)
	    os << "]:" << it.value();
	    //it.index(); // inner index, here it is equal to it.row()
		os << ", " ;
	  }
	os << " }";
	return os;
}

namespace LuaNumeric {

template<class T>
sol::object getter(const Eigen::SparseMatrix<T>& vec, sol::table key, sol::this_state L) {
		// we use stack_object for the arguments because we know
		// the values from Lua will remain on Lua's stack,
		// so long we we don't mess with it
		// String keys failed, check for numbers
	sol::object key1 = key[1];
	sol::object key2 = key[2];

	auto maybe_numeric_key1 = key1.as<sol::optional<int>>();
	auto maybe_numeric_key2 = key2.as<sol::optional<int>>();
		if (maybe_numeric_key1&&maybe_numeric_key2) {
			int n1 = *maybe_numeric_key1;
			int n2 = *maybe_numeric_key2;
			if( n1 <= vec.rows() && n1 >= 1
				&& n2 <= vec.cols() && n2 >= 1)
				return sol::object(L, sol::in_place, vec.coeff(n1-1, n2-1));
		}
		// No valid key: push nil
		// Note that the example above is a bit unnecessary:
		// binding the variables x and y to the usertype
		// would work just as well and we would not
		// need to look it up here,
		// but this is to show it works regardless
		return sol::object(L, sol::in_place, sol::lua_nil);
}
template<class T,class T1>
void setter(Eigen::SparseMatrix<T>& vec, sol::table key, const T1& val) {
		// we use stack_object for the arguments because we know
		// the values from Lua will remain on Lua's stack,
		// so long we we don't mess with it
		// String keys failed, check for numbers
	sol::object key1 = key[1];
	sol::object key2 = key[2];

	auto maybe_numeric_key1 = key1.as<sol::optional<int>>();
	auto maybe_numeric_key2 = key2.as<sol::optional<int>>();
	if (maybe_numeric_key1&&maybe_numeric_key2) {
		int n1 = *maybe_numeric_key1;
		int n2 = *maybe_numeric_key2;
		if( n1 <= vec.rows() && n1 >= 1
			&& n2 <= vec.cols() && n2 >= 1){
			vec.coeffRef(n1-1, n2-1) = val;
			return;
		}
	}
	throw "Eigen::SparseMatrix out of range or index error";
}
	typedef std::complex<double> complex;
	void open_luaSparseEigen(sol::table& module) {
		module.set_function("new_sparse_matrix", sol::overload(
				[]() { return SpMat(); }
				,[](int m, int n) { return SpMat(m,n); }
				)
		);
		module.set_function("new_sparse_cmatrix", sol::overload(
						[]() { return SpCMat(); }
						,[](int m, int n) { return SpCMat(m,n); }
						)
				);

		module.set_function("eigs", sol::overload(
					[](SpMat& A,int nvec) {
						Eigen::VectorXd v;
						int res = solve_sp_eig(A,v,nvec);
					},
					[](SpMat& A,SpMat& B,int nvec) {
						Eigen::VectorXcd lambda;
						Eigen::MatrixXcd V;
						std::tuple<int,Eigen::VectorXcd,Eigen::MatrixXcd> res;
						std::get<0>(res) = solve_sp_geig(A,
													  B,
													  std::get<1>(res),
													  std::get<2>(res),
													  nvec);
						return std::move(res);
					}
/*

					,[](SpCMat& A,SpCMat& B,int nvec) {
						Eigen::VectorXcd lambda;
						Eigen::MatrixXcd V;
						std::tuple<int,Eigen::VectorXcd,Eigen::MatrixXcd> res;
						std::get<0>(res) = solve_geig(A,
													  B,
													  std::get<1>(res),
													  std::get<2>(res),
													  nvec);
						return std::move(res);
					}
//*///////////
				)
		);

		module.set_function("eigs_spectra", sol::overload(

					[](SpMat& A,SpMat& B,int nvec) {
						Eigen::VectorXcd lambda;
						Eigen::MatrixXcd V;
						std::tuple<int,Eigen::VectorXcd,Eigen::MatrixXcd> res;
						std::get<0>(res) = solve_gsheig_spectra(A,
															  B,
															  std::get<1>(res),
															  std::get<2>(res),
															  nvec);
						return std::move(res);
					}
				)
		);

		module.new_usertype<SpMat>("SpMat"
				,sol::no_constructor
				//sol::constructors<EigenMatrix(EigenMatrix::Index, EigenMatrix::Index)>(),
				,"resize", (void(SpMat::*)(SpMat::Index, SpMat::Index)) &SpMat::resize
				,"get", [](SpMat* m, int nm,int nn )->double& {
							return (*m).coeffRef(nm-1, nn-1);
				}
				,"set", [](SpMat* m, int nm,int nn, double v){
						(*m).coeffRef(nm-1, nn-1) = v;
				}
				,"set_row", [](SpMat& m, int nm,double v){
					if(nm<1 || nm > m.rows()){
						throw " out of range ";
					}
					for( int ii=0; ii < m.cols(); ii++){
						m.coeffRef(nm-1, ii) = v;
					}
				}
				,"set_col", [](SpMat& m, int nm,double v){
					if(nm<1 || nm > m.cols()){
						throw " out of range ";
					}
					for( int ii=0; ii < m.rows(); ii++){
						m.coeffRef(ii,nm-1) = v;
					}
				}
				,"set_identity",[](SpMat& m){
						m.setIdentity();
				}
				,"set_zero",[](SpMat& m){
						m.setZero();
				}
				,"solve", [](SpMat* m, const Eigen::VectorXd& b, Eigen::VectorXd& x){
					int res = solve(*m,b,x);
					return res;
				}
				,"solve_lu", [](SpMat* m, const Eigen::VectorXd& b, Eigen::VectorXd& x){
					int res = solve_sp_lu(*m,b,x);
					return res;
				}
				,"solve_ldlt", [](SpMat* m, const Eigen::VectorXd& b, Eigen::VectorXd& x){
					int res = solve_sp_ldlt(*m,b,x);
					return res;
				}
				,"solve_cg", [](SpMat* m, const Eigen::VectorXd& b, Eigen::VectorXd& x){
					int res = solve_sp_cg(*m,b,x);
					return res;
				}
				,"solve_eig", sol::overload(
						[](SpMat& M, Eigen::VectorXd& v,int nvec){
							int res = solve_sp_eig(M,v,nvec);
							return res;
						}
						,[](SpMat& M, Eigen::VectorXcd& v,int nvec){
							int res = solve_sp_eig(M,v,nvec);
							return res;
						}
				)
				,"prune", [](SpMat& M)
				{
					M.prune(0.0);
				}
				,"rows", &SpMat::rows
				,"cols", &SpMat::cols
				,sol::meta_function::index, &getter<double>
				,sol::meta_function::new_index, &setter<double,double>

		);

		module.new_usertype<SpCMat>("SpCMat"
				,sol::no_constructor
				//sol::constructors<EigenMatrix(EigenMatrix::Index, EigenMatrix::Index)>(),
				,"resize", (void(SpCMat::*)(SpCMat::Index, SpCMat::Index)) &SpCMat::resize
				,"get", [](SpCMat* m, int nm,int nn )->std::complex<double>& {
							return (*m).coeffRef(nm-1, nn-1);
				}
				,"set", sol::overload(
						[](SpCMat& m, int nm,int nn, std::complex<double> v){
						m.coeffRef(nm-1, nn-1) = v;}
						,[](SpCMat& m, int nm,int nn, double v){
							m.coeffRef(nm-1, nn-1) = v;}
				)

				,"set_row", sol::overload(
					[](SpCMat& m, int nm,double v){
						if(nm<1 || nm > m.rows()){
							throw " out of range ";
						}
						for( int ii=0; ii < m.cols(); ii++){
							m.coeffRef(nm-1, ii) = v;
						}
					},
					[](SpCMat& m, int nm, const std::complex<double>& v){
							if(nm<1 || nm > m.rows()){
								throw " out of range ";
							}
							for( int ii=0; ii < m.cols(); ii++){
								m.coeffRef(nm-1, ii) = v;
							}
					}
				)
				,"set_col",sol::overload(
					[](SpCMat& m, int nm,double v){
						if(nm<1 || nm > m.cols()){
							throw " out of range ";
						}
						for( int ii=0; ii < m.rows(); ii++){
							m.coeffRef(ii,nm-1) = v;
						}
					},
					[](SpCMat& m, int nm,const std::complex<double>& v){
						if(nm<1 || nm > m.cols()){
							throw " out of range ";
						}
						for( int ii=0; ii < m.rows(); ii++){
							m.coeffRef(ii,nm-1) = v;
						}
					}
				)

				,"prune", [](SpCMat& M)
				{
					M.prune(std::complex<double>(0.0));
				}
				,"set_identity",[](SpCMat& m){
						m.setIdentity();
				}
				,"set_zero",[](SpCMat& m){
						m.setZero();
				}
/*
				,"solve_eig", sol::overload(
						[](SpCMat& M, Eigen::VectorXcd& v,int nvec){
							//int res = solve_eig(M,v);
							//return res;
						}
				)
*/
				,"solve", [](SpCMat* m, const Eigen::VectorXcd& b, Eigen::VectorXcd& x){
					int res = solve(*m,b,x);
					return res;
				}
				,"solve_lu", sol::overload(
						[](SpCMat& M, const Eigen::VectorXcd& b, Eigen::VectorXcd& x){
							int res = solve_sp_lu(M,b,x);
								return res;
						}
						,[](SpCMat& M, const SpCMat& B, SpCMat& X){
							int res = solve_sp_lu(M,B,X);
							return res;
						}
				)
				,"solve_ldlt", [](SpCMat* m, const Eigen::VectorXcd& b, Eigen::VectorXcd& x){
					int res = solve_sp_ldlt(*m,b,x);
					return res;
				}
				,"solve_cg", [](SpCMat* m, const Eigen::VectorXcd& b, Eigen::VectorXcd& x){
					int res = solve_sp_cg(*m,b,x);
					return res;
				}

				,"rows", &SpCMat::rows
				,"cols", &SpCMat::cols
				,sol::meta_function::index, &getter<std::complex<double>>
				,sol::meta_function::new_index, sol::overload(
						&setter<std::complex<double>,std::complex<double>>
						,&setter<std::complex<double>,double>
				)
				,sol::meta_function::multiplication, sol::overload(
							[](const SpCMat& a, const SpCMat& b){
								SpCMat c = a*b;
								return std::move(c);
							}
				)
				,sol::meta_function::addition, sol::overload(
							[](const SpCMat& a, const SpCMat& b){
								SpCMat c = a+b;
								return std::move(c);
							}
				)
				,sol::meta_function::subtraction, sol::overload(
							[](const SpCMat& a, const SpCMat& b){
								SpCMat c = a-b;
								return std::move(c);
							}
				)

		);

	}
} // namespace my_object




